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Abstract 


A linear analysis and comparison of the damping properties of six dynamic 
Initialization schemes Is presented, Indicating that the Okamura-Rlvas scheme 
has the most efficient damping properties over the whole frequency range, and 
suggesting that it should be faster than the other methods and give more stable 
results. The results obtained with a nonlinear shallow water equations model 
agree well with the linear analysis. The Okamura-Rlvas scheme attains complete 
balance In the equivalent of 3 to 6 hours of leapfrog forecasting, and requires 
In this model an order of magnitude less computation than the balance equation 


solution. 
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1. Introduction 

It 1b well known that the use of observed data directly as Initial fields 
for numerical primitive equation models results In predicted gravity-inertia 
oscillations with amplitudes far in excess of those observed In the atmosphere. 

These oscillations may mask the meteorologically significant motions, as was 
first observed In Richardson’s (1922) experiment. The problem Is that the ob- 
served fields contain excessive Imbalances between Coriolis and pressure forces, 
which arise primarily from observational errors. 

Two techniques have been widely used to solve the Initialization problem, 
l.e., to redufv the excessive Imbalances 'between the mass and wind fields. 

The fivst :^8 the balance equation approach, In which the initial mass and 
wind fields are forced to satisfy some form of the balance equation, from a 
simple geostrophlc approximation to the full nonlinear balance equation. Phillips 
(1960) suggested that a consistent estimate of the divergent wind should also be ob- 
tained from the quasl-geostrophlc omega-equation and added to the nondlvergent 
wind determined through the balance equation. H^aghton and Washington (1969) 
and Houghton, Baumhefne:: and Kasahara (1971) show how this method can be 
applied on a global scale. 

This method is generally satisfactory, but It has disadvantages, some 

of which are: a) It requires the numerical solution of a set of "diagnostic" 

type of equations which may not be directly compatible with the prediction 
equations; b) In extratroplcal regions, where ;;lnds are usually obtained In 
terms of the mass field, the nonlinear balane** equation is of a mixed hyperbolic- 
elliptic type, and Is readily soluble only when purely elliptic. Non-elliptlcity 
occurs in strongly anticyclonic regions where the total vorticity Is smaller 
than one half of the Coriolis parameter. This problem, which as we will see 
becomes more critical the better the resolution of the model, is usually "solved" 
by artificially modifying the observed mass fields in non-elliptical regions, 
but obviously this is not a satisfactory procedure. 
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The second technique which Is becoming increasingly popular is known as 
"dynamic Initialization". In this technique the primitive equations themselves 
are used to march forward one or several time steps and then return to the Initial 
time by reversing the time step. A damping scheme, as for example Euler' back- 
wards, is used, and If the procedure is repeated, the high-frequency Inertia - 
gravity waves generated by the initial Imbalance are eventually damped out 
(Nltta and Hovermale, 1969) . Nltta and Hovermale suggested that the original 
mass field be recovered after each forward-backward iteration under the assump- 
tion that, at least in extratropical regions, measurements of pressure and height 
fields are more reliable than winds. Partial recovery of the height fields has 
been proposed by Mesinger (1972) and Winninghoff (1973). 

The advantages of the dynamic initialization techniques are two: a) 

they are very simple to use. The forecast equations (minus the dissipative terms) 
are used, and therefore the initialized fields are compatible with the forecasting 
scheme, b) there is no elllpticlty constraint when winds are obtained from the 
mass field. On the other hand, they have a serious disadvantage which is their 
relative inefficiency. The original Nitta-Hovermale scheme may require several 
hundred iterations before convergence, equivalent to the number of computations 
necessary to make a few days' forecast, 

A similar, but somewhat more efficient scheme, has been proposed by Mesinger 
(1972). Temperton (1973) has proposed an initialization scheme in which the 
fields obtained after forecasting 6 timesteps, say, starting from the initial 
time, and those obtained by "hlndcasting" 6 time-steps starting again from the 
initial time, are averaged, and the procedure is then repeated. The method is 
based not on the damping properties of the time scheme, but on the assumption that 
while large scale meteorologically significant components of the fields will 
remain largely unchanged, gravity waves generated by the imbalances will be 
averaged out. 


Okamura (see Nitta, 1969, appendix) has proposed a scheme of the Nitta- 
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Hovemale type, which is simpler and more efficient, but has not found wide 
use. Kalnay de Kivas has developed a more flexible version of Okamura's scheme 
which enhances its effectivity. Increases its stability and still retains its 
simplicity. 

The purpose of this paper is to compare the balance equation method and 
the different d^’namlc initialization techniques. Section 2 contains a brief 
discussion of the solution of the balance equation. In section 3 we discuss the 
mathematical properties of the different dynamic initialization methods. The 
design of the numerical experiments is presented in section 4 and section 5 
contains the numerical results. 
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2. The balance equation approach 

The balance equation, as proposed by Chamey (1955), Is obtained from the 
divergence t juatlon by neglecting divergence and vertical motion terms, and by 
assuming a non-dive rgent wind, v » k x Vd;. This yields 

V^(|. - -Vtft-vf + 2\l>^ - 2\ii ^ •= 0. (2.1) 

xy XX yy 

This equation ’has a maximum relative error O(RlRo), where R1 and Ro are 
the Richardson and Rossby numbers, and therefore Is accurate at least up to 
10% even In equatorial regions. 

If the stream function is considered known from the observed wind field, 
the balance equation can be solved to obtain the geopotential and Is obviously » 

elliptic. On the other hand, If <f> Is taken directly from the data, and 
(2.1) Is solved for i|i, It Is easy to show that (2.1) Is elliptic If 


Vf • f 

f " f " 2 ' 


( 2 . 2 ) 


It should be emphasized that this restriction is non-physical, and does 
not arise when the complete consistent system of balance equations is consi- 
dered (Charney, 1973) . Houghton and Washington (1973) have shown that the 
solution of (2.1) for <|) yields more accurate results In the tropics, wh:ll.i the 
solution for i{) is preferred for extratroplcal regions. Thus the elllpticity 
condition may be violated In isolated regions over a substantial portion of the 
globe. 

In extratroplcal regions the balance equation is usually solved by rewriting 
eq. (2.1) after Betterssen (1953) and Bolin (1955) as follows; 

V^tp = -f ± (2V^(f> + f2 + _ 2vf . Vip)^ (2.3) 


where A = and B = ^'(^xy deformation terms. The positive or 

negative sign of the radical is applicable in the northern or southern hemisphere 
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respectively. A procedure for solving eq. (2.3) for ip, called the cycle-scan 
method by Mlyakoda (1956 » 1960) and Schuman (1957) couslsts of evaluating the 
right hand from a given guess for and correcting f by Inverting the 

laplaclan on the left hand side. The pro e usually converges If the elllptlclty 
constraint is satisfied, which also ensures that the right side of 2.3 Is real. 


« 
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3. Analysis of the dynamic Initialisation methods 

In this section we make a linear analysis and comparison of the properties 
uf six dynamic Initialisation methods. For simplicity, let us represent a 
set of primitive equations by 


la -Fu (3.1) 

O t 

where u is a column vector of the dependent variables and^ is a matrix differen~ 
tial operator which contains no explicit time derivatives. This equation can be 
converted into Its flni'^e difference (or ’spectral) equivalent. For example, an 
Euler marching scheme may be represented as ^ 

^(t+I) ^ ^(t) ^ « (I + AtF)U^'^^ (3.2) 


where arc now discrete values of the dependent variables at time t = xAt 

and F is the finite difference (or spectral) equivalent ofT”. 

The first two schemes we will discuss were suggested by Nltta and Hovermale 
(1969). The first method, which we will call Nitta-Hoverraale 1 or NHl, is the 
following: 


= U 

u = u 

. u 


+ AtFU 
+ AtFU 

- AtFU 

- AtFU 


(v) 

* 

(p) 

AA 


(3.3a) 


(3.3b) 


Note that it consists of a single forward time step using the Euler back- 
ward or Matsuno scheme (3.3a) followed by a backward time step using the same 
scheme (3.3b). The superindex (v) indicates number of iterations, not time, 
since at the end of a complete Iteration we are at the same initial time level. 

A recursion relation involving succesive values of U is obtained from 


(3.3): 
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y(*+l) , ^ ^ At2p2 + 


(3. A) 


The second method (Nltca-Hovermale 2 or NH2) uses a modified Euler-backward 


scheme: 


U 


U 

u 

(y) 

u 


u 


u 

(vfl) 


- + AtFU* 

- + titFU** 

. u(w) _ |t pu<y) 

- - AtFU , 

= - AtFU. 


(3.5a) 


(3.5b) 


The recursion relation for succeslve Iterations is 


u(v+i) . .J ^ ^ p6) 


(3.6) 


A third method was proposed by Okamura (see Nitta, 1969, Appendix). It 
Involves an explicit Euler time scheme applied forward and backward , 


U* = + AtFU^'’^ 


it* * * 

U = U - AtFU 


(3.7a) 

(3.7b) 


✓ X 

followed by a linear combination of U'^'' and U : 


. 3U<'’> - 2u“ 


(3.7c) 


A fourth method is a generalization of Okamura 's method by Kalnay de Rivas, 
in which (3.7c) is replaced by 


!)<'*» = (n + 1) U<'’> - „u‘* 


(3.7d) 


and n is allowed to vary with v . It reduces to Okamura *s scheme if n = 2, 
The recursion relation for the Okamura-Rivas (0-R) scheme is 
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u(v+l) . (I +nAt2p2) 


(3.8) 


Before Introducing Che last two methods, let us examine the stability and 
damping characteristics of these four schemes. For this purpose, consider a 
single harmonic wave of frequency tu In time, Iji « U e^^. u; represents any of 
the characteristic frequencies of the problem (Rossby waves, inertia-gravity 
waves, etc.), or, more precisely, the eigenvalues of the linearized version of 
the operator F, and U are the corresponding eigenvectors. Then, Che time dif- 
ferencing can be expressed explicitly: 


AtFU « ImAttl « ipU. 


(3.9) 


Introducing (3.9) into 


(3.4) 

U 


, (3.6) and (3.8), we obtain for each method 




(3.10) 


where R is the damping factor corresponding to one complete iteration and Is 
given by 



R •= 1 - p2 + p4 

NHl 



R o 1 - p2 + Jjjpfi 

NH2 

(3.11) 


R = 1 - 2p^ 

0 



R = 1 - np2 

0-R. 


Stability of 

the iterative methods requires 

|r| <1 and this places 

a 

restriction on the 

size of the time increment At: 




At^ < l/o)2 
— m 

NHl 



At^ < 2/(i)2 

NH2 

(3.12) 


- m 



At2 < 2/(ntij2) 

— m 

0-R . 



Here is ! ■naximum frequency present in the problem (for example the 
frequency of the shortest Lamb wave in a primitive equations model) . We have 
assumed in (3.12) that n is fixed, and in that case it is clear that Okamura’s 
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choice of n 2 maximizes the efficiency of the method. 

Fig. 1 shows the absolute value of the damping factors |R| as functions 
of p for one complete Iteratlony as In (3.11). The damping of one 0-R Iteration 
Is shown for n * 1, 2 and 4. Note that NHl and 0-R for n ■> 2 (Okamura's scheme) 
show the undesirable property that the highest frequencies (p 'v 1) are not 
damped, and 0-R for n ■ 4 is obviously unstable, even though it is strongly 
damping for low and mid-range frequencies. 

In the 0-R method, this can be corrected by allowing n to repeatedly take 
on a sequence of values during the Iterative process. The total damping factor 
of a sequence is the product of the damping factors at each n of the sequence. 

The sequence can be chosen in such a way as to maximize the damping both at mid- 
range and high frequencies. A good example of such a sequence, although not 
necessarily an optimum one, is n = 1, 1.6, 4. The numbers in this sequence were 
specifically chosen to correct for the sharp rises in the single-n curves shown 
in Fig. 1, 

It should be noted that the number of computations involved in a single 
complete iteration depends on the method used. The most costly part of the 
scheme is the computation of each operation F on U. A meaningful comparison of 
the cost in computer time required by the methods is therefore the number of 
F computations required per complete iteration. Another factor that may deter- 
mine the advantage of a method is the amount of computer memory required, given 
approximately by the number of sets of the dependent variable U which must be 
stored separately. Table 1 comj)ares the four methods, as well as Mesinger's 
and Teraperton's methods, on the basis of these practical considerations. 

Fig. 2 shows the damping of the NHl, NH2, Okamura and Okamura-Rivas with 
n = 1, 1.6, 4, after 12 F computations, and therefore compares the relative 
efficiency of the methods. Note that the 0 and 0-R methods are far more efficient 
than the NHl and NH2 methods for low and mid-range frequencies, but only the 0-R 
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Method 

if of F per complete iteration 

y/U stored 

NHl 

4 

2 

NH2 

6 

3 

Okamura 

2 

3 

Okamura^Rivas 

■ 2 

3 

Mesinger 

4 

2 

Temperton (N) 

2N 

3 


Table 1; Comparison of number of computations and storage requirements of 


the iterative schemes. 
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method exhibits strong damping of high frequencies. 

Meslngcr (1972) has proposed a marching scheme which is a generalization 
of the Matsuno and Heun schemes: 

U* . -f aAtFU^'’^^ 


u(T+l) 


U<t) 


* 

+ AtFU . 


(3.13) 


The extrapolation factor ^ ■ 1 corresponds to the Matsuno scheme, and ^ 
coincides with the Heun scheme. Meslnger suggests the use of thl>' scheme in a 
forward -backward fashion for dynamic initialization. After a complete iteration, 
the damping factor Is 


R - 1 - (2a - l)p2 + p**. (3.14) 

The time step At is chosen by Meslnger as 4t = (2a - l)/(2a^o}^) so that 
damping Increases ir.'tiiiiLonically with frequency. Replacing At^ in (3.14) we 
obtain 

R - 1 - 2(1 - %a) (^)^ + (1 - i£a)^ (^)*^ . (3.15) 

m m 

Then maximum damping is obtained by letting 


- ^ v2 ,tu 

®MES “ ^ ” 2(~) + (— 


(3.16) 


m m 

Note that this method , like NHl requires 4 F computations to complete one 
iteration. 

A distinct advantage of Meslnger 's scheme is that it can also be used as 
a regular forecasting scheme, in the same way as Matsuno's scheme. However, 
in order to provide significantly more damping than Matsuno’s scheme, a^ has 
tc be much larger than one. The time truncation error of this scheme is 


T 


(h -a) At 


3^u , At^ 3 ^u 
^ 6 3t^ 


(3.17) 


i.e., first order in At (unless ^ * ^) and proportional to Meslnger suggests 
the use of aAt = 6 hours, but we think this will introduce intolerably large 
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errors even for sy '' optic waves. 

Fig > 3 compares the damping of the Mesinger and Okamura->Rlvas schemes 

after 12 computations of F. In the 0-R scheme, At - l/«^* so that p * u)At ■ ^ . 

m 

Even though Meslnger's scheme Is better than the NHl and NH2 methods, It Is 
considerably less efficient <than the 0-R scheme with n ■■ 1, 1.6, 4. 

The last scheme that we will discuss is the one proposed by Temperton 
(1973) . He suggests that a forward Euler time step followed by centered leap 
frog steps be applied so as to obtain . Then, starting again from the 

initial time a similar forecast with negative time step should be executed, 
obtaining An iteration is completed by averaging the two forecasts: 

^(v+1) ^ j^jy(NAt) _^*y(-NAt)^ (3.18) 


and the process can then be repeated. To understand how this method works, 
consider the wave equation 3u/9t + c3u/3x« 0, and suppose that initially « consists 
of a single "bump” of width If we apply Temperton *s method once, with 
T ■» NAt, assuming that we maki. no truncation error, we obtain 


u^^^ = Jslu*^^^(x - cT, 0) + u^^^(x + cT; 0)]. 

If cT << L, u^^^* u^^\ If cT > L, u^^' will contain two bumps, with one 

half of the amplitude of the initial bump, and therefore the energy contained 

in u is now halved. This analysis suggests that T should be chosen in such a 

way that c T « L for synoptic waves , and c T > L for gravity waves . The 
S3 g ~ 8 

optimum T found experimentally by Temperton, T = 6At, satisfies this criterion 

for the values of c “ 140 m sec ^ and L = 2Ay = 400 km corresponding to his 

g g 

numerical model. 

If we apply Temperton ’s method to a single harmonic component U 
with N = 6, we obtain a total damping 


Rg = 1 - 18p^ + 48?** “ 32p®. 


(3.19) 
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This Involves 12 computations of F and is also plotted on Fig. 3. Note 

> that If At is taken as l/u> , the maximum value allowable for stability » as done 

m 

in Fig. 3, then the damping for low frequencies of Temperton's scheme is 

i 

' Rg 1 - 18p^. The 0-R scheme damping after 12 F computations is * 1 - 13. 2p^. 

However, from "Fig. 3 we see that Temperton's scheme does not provide damping 
at all frequencies unless the time step is taken as At < 1/ (2o)^) . In that case 
Temperton's method becomes much less efficient than the O-R method e\*ca at low 
frequencies. It should also be noted that for odd N, Temperton's scheme will 
amplify high frequencies. 

In summary we have shown from a linear analysis that of the six methods 
compared, the Okamura-Rlvas scheme with n » 1, 1.6, 4, which is not necessarily 
optimum, provides the most efficient damping at low and middle frequencies, and 
the best damping at high frequencies. 

Before closing the section we want to make two comments: 

a) We have assumed that (o is a real number. This implies that irreversible 

terms, like friction, are not Included, a..d amplifying or decaying modes, if 

2 

present, need to have slow growth rates . 

b) In this analysis we have not allowed for the restoration of heights 
after each Iteration. 


2 

Kalnay de Rivas (1975) has recently developed a simple scheme that will damp 
out any time dependent behaviour, and can be used, for example, to obtain an 
unstable steady state solution of the complete primitive equations, including 
forcing and dissipative terms. 



4# Design of the numerical experiments 


To test the various initialization techniques, a nonlinear, shallow 
water model was developed, defined by the following equations; 


3 ^u 

at 


- iiH. + f - 4, 

3 x 3 y 'T r 


dx 


(4.1a) 


34>v 

9t 


liiLV _ ^ li 

3x 3y ^ ^ 3y 


(4.1b) 


91 

3t 


34iu _ 3(j)V 

3x ” 3y 


(4.1c) 


where (u,v) are the horizontal velocity components, 41 = gh, and h is the height 
of the free surface. The Co^olls parameter is held constant, f >* 10 '♦s 

In order to avoid nonlinear instability, avi energy conserving finite dif- 
ference scheme is used for the right hand side of (4.1). The grid used has 
4x = Ay ■ As = 250 km. 

In our experiments we artificially generate an initial state which is 
as closely in balance as possible. Then the velocity and height fields of the 
balanced data are altered or perturbed by various methods, and the different 
initialization techniques used to restore bal-vnce. We are interested in deter- 
mining how well the initial "correct" fields are recovered, and how well in 
balance the initialized fields are. 

The initial balanced state is generated by integrating the model for a finite 
time T with an artificial mass source term S(x,y,t) added to equation 4.1c. The 
integration is started from a state of rest with a level free surface at H = 3 km. 

A leapfrog marching procedure is used with a forward . ?ler step every 24 leapfrog 
steps to avoid the separation of fields at even and odd steps. The spatial 
variation of the source function is 

S(x,y,t) = S(t) sin(|^ x) sin(^ y) 
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where L > 4000 km. It Is convenient to visualize the state produced as a checker- 
board pattern extending periodically over the Infinite £-plane. Numerical compu- 
tations are nade on a rectangular grid which ext nds one wavelength In the x 
direction and half a wavelength In the y direction. Boundary conditions are 
periodic, with the north-south boundaries matched diagonally to preserve periodicity. 

We found that the degree of balance attained by the Initial state was 
critically dependent on the time variation of the source term, S(t). The field 
shown In Fig. . 4 Is attained with the following form 

S(t) - 1“ sln(~ t) 

where the Integration was terminated at time T = 8 days, and A ■ 1.01 x loV^s. • 

Is the Integrated strength. The contours are drawn at Intervals of 60 m; the 
low Is 340 m below the mean height, and the high 150 m above It. The associated 
velocity fields (not shown) are cyclonic around the low and antlcyclonlc around 
the high with 30 m s - ^ maximum speed. 

It was found necessary to add the source In very small Increments (At = 5 min) 

In order to minimize the generation of Imbalances In the Initial state. After 
the source was "turned off’, the forecast was continued using again the leap- 
frog scheme with At = 12 min. When the height at the point P (see Fig. 4) was 
monitored, it was found that the amplitude of inertial-gravity waves present was 
only 0.2 m, indicating a high degree of balance. Similar attempts to generate 
Initial data were made using an exponential and a linear time dependence for 
S(t), with the same Integrated strength. The resulting forecasts, after the 

source was turned off, showed gravity waves with amplitudes 25 and 100 times 

3 

greater , respectively . 

^hls observation can have an important application in numerical climate studies, 
in which a physical parameter P (sea surface temperature, ground albedo, etc) 
is changed from its standard value Pq to a different value P. , to study the 
effect that this change has on the circulation. Our observations suggest that 
the initial Imbalance, and consequently the amplitude of the transient inertia- 
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The Initial fields (fig. 4) and the fields obtained after a 48 hours £ore<* 
cast are used as standards In computing rms departures of the velocity and 
height fields obtained using the Initialization methods. In the Initialization 
methods we allowed the option of restoring the geopotential heights after eanh 
iteration. 


« 


(eont) gravity waves generated by the change, are maximized by the common prac- 
tice of changing Pq to abruptly. On the other hand, a smooth, change over a 
time T of the form 


P(t) 


pq +n 
2 


, Pq - Pi 
2 


cos 


TTt 

T 


0 < t < T 


will minin^ize the generation of gravity waves. When P is changed abruptly, P(t) 
is a Heaviside function, and its Fourier transform contains large amplitude 
high frequencies, which therefore excite gravity waves. When P is changed 
smoothly as we suggest, P and P’ are continuous for 0 < t < T, and if T is of 
the order of a few days, the transform of P(t) will contain only small ampli- 
tude high frequencies. We acknowledge a useful discussion with Mr. Mark Cane 
about this point. 
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5. Numerical results 

As Indicated in the previous section, a balanced "synoptic wave" with a 
geopotential amplitude of 250 m and maximum winds of 30 m s ^ was obtained 
and used as reference state for the Initialization experiments. 

In order to Introduce Initial imbalances representing the effect of ob- 
servational errors, the reference fields were perturbed in two ways. First we 
assumed that the "observed" geopotential field contained no errors, but replaced 
the velocity field by the geostrophlc velocity field. Second, we simulated 
observational errors by adding to the balanced fields normally distributed ran- 
dom nu'jibers , with standard deviations typical of atmospheric measurement errors. 

The damping of each Iterative met^od was* maximized by using the maximum 
time step (in whole minutes) for which the method remained stable. TaOi-e 2 
compares these experimental values to the upper limit from the linear stability 
criteria as expressed in (3.12). 

Hote chat the iterative schemes adhere to the linear criteria more closely 
than does the leapfrog scheme, and that for the NH2 method, there Is strong 
damping of high frequencies even though ta At > 1 (fig. 2). 

a) Initialization of the geostrophically perturbed state 

The geostrophically determined winds depart from the reference wind fields 
by an ms error of 7.7 m s If the forecast proceeds without initialization, 

gravity waves with amplitude of 125 m appear, which are sufficient to strongly 
distort the synoptic wave which has an amplitude of 250 m. 

If the ellipticity condition (2.2) is satisfied everywhere, the balance 
equation can be solved without altering the geopotential field. Our reference 
state has enough amplitude that (2.2) is violated at a few points around the 
perimeter of the high. At the points where 


ij 


, 4Sj4-ii4i + fi 


As 


< 0 . 


(5.1) 
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Method 

At (nonlinear) 

At (linear) 

NHl 

16 min 

17 min 

NH2 

22 min 

, 24 min 

0 

16 min 

17 min 

0-R 

17 min 


Leapfrog 
(forecast and 

12 min 

17 min 

Temperton’s 

method) 

•* 


Table 2: Maximum time 

steps according to the linear 

criteria and in the 


nonlinear numerical model 





we applied the following "correction". We defined =• from 

the conditions A + B =» 1, 


7 ,/. u ) + ii = 

As^ 2 



we determined A and B, Here, 1 ( 1 ^^ is the average geopotential at the four 
adjacent points, and a is a small non->negative number. After replacing i|i 


ij 


I 


by (|i^j at the points where < 0 we repeated the process until > 0 

at all points. In our case, the best results were obtained with a 0, requiring 

5 iterations. The height field was modified at 8 points, resulting in a maxi- 

4 

mum change of 0.5 m. 

Once the geopotential satisfied the ellipticity condition everywhere, 

28 cycle-*3can iterations were sufficient to solve the balance equation. Table 
3 shows that the application of the balance equation reduces the rms wind error 
to 0.7 m sec and the amplitude of the residual gravity wave to 3m. We conclude 
that, in this case, the balance equation does provide a significant restoration 
of both the initial fields and their state of balance. 

In additional experiments, however, we found that when the ellipticity 
condition was violated at a larger number of points the ^‘correction” procedure 
became insufficient to “elliptize” the geopotential field. When the source 
strength was increased by 20% the correction procedure did not converge, and 
therefore we were unable to solve the balance equation . The iterative methods, 
of course, were not affected by the change in source strength. We will come 
back to this point in subsection 5c. 


^ An anonymous reviewer has pointed out the existence of a better “correction” 
procedure: Calculate h = V^iJ) + f^/2 at every point. Then solve for 

+ f^/2 = h^, where h- r h if h > 0, h^ If h- 0, and e. is a small positive 
number. will satisfy the ellipticity condition. 


In the first series of dynamic initialization experiments we performed 150 
complete iterations with the NHl, NH2, 0 and 0-R schemes, restoring the heights 
after each iteration. If we assume that each evaluation of the time derivative F 
(eq. 3^2) is equivalent to one timestep in a forv;ard forecast, the equivalent 
time traversed in 150 iterations is 5 days for the NHl method. 
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7.5 days for the NH2 method » and 2.5 days for the 0 and 0-R methods. 

Fig. 5a shows the decrease In rms error with the number of Iterations. 

A considerable reduction of error occurs, but even 150 iterations are not suffi- 
cient for any method to converge to a steady value of the rms error. Figures 
6a, b, c and d (lower curves) compare the forecasts of the height at the point 
P, when the heights are restored during initialization. Clearly, even 150 
Iterations are not sufficient to restore balance when the heights are not allowed 
to adjust freely. This can be understood in the following way: after the first 

iteration, the heights are modified from the original values with an rms change 
of 10-20 m, depending, on the method. When we restore the geopotential field 
to correct this change, we are restoring a large portion of the initial imbalance . 
By restoring the heights after each Iteration we force the slow convergence of 
the iterative methods. 

In the next experiments we performed 150 iterations allowing the heights 
to adjust freely. Figures 5b and 5c show the variation of the zms wind and 
height errors with the number of iterations. Only the first 30 iterations 
are shown. All four methods reached a steady rms wind error of 6.9 m s ^ 
and created an rms height error of 46 m. To compare the convergence rates, 
we estimate the number of iterations required to reach the steady values: 

40 (NHl) , 15 (NH2) , 15 (O) and 12 (O-R) . The equivalent time traversed during 
these iterations is 32 h (NHl) , 18 h (NH2) , 6 h> (0) and 5 h 
(0-R). Figures 6a, b, c and d (upper curves) show the height forecast at the 
point P when the heights adjusted freely during the 150 iterations. They show 
that high frequency oscillations with increasing amplitude appear in the forecasts 
after using NUl and 0. These are due to the aliasing of very high frequency 
waves not eliminated by the NHl and 0 methods. On the other hand the NH2 and 
0-R methods completely eliminate these oscillations if the height is allowed 
to adjust. As shown in fig. 7, only 12 to 15 0-R iterations, or the equivalent 
of 5 to 6 hours of forecast, are enough to attain balance and eliminate inertia- 




ID' 


Table 3: Initialization of the geostrophically perturbed case 
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gravity waves from the forecast. 

We made several experiments with Temperton's method. In the first we 
performed 150 complete Iterations with N 6, l.e. marched forwarrf from the 
Initial time 6At» marched backward from the Initial time 6 At, averaged the winds 
and restored the Initial heights. Table 3 shows that the reduction of error Is 
smaller than for the 0 and 0-R methods, but there Is more balance, as shown by 
the smaller amplitude of the residual Inertla-gravlty waves. Xt should be noted, 
however, that this method Is 6 times more expensive In computer time than the 
0>'R method. In fig. 5a the rms error Is shown by the crosses plotted every 6 

I 

0-R Iterations. Clearly the 0-R method is much more efficient than Temperton’s 
method. • 

When we allowed the heights to adjust freely during 30 Temperton iterations 
(equivalent to 180 0-R iterations), the rms errors were similar to those of the other 
schemes, but gravity waves with amplitudes of 2-15 ra were present in the forecast; this 
agrees with our analysis in section 3, in that Temperton's scheme does not filter 
out high frequency waves uniformly. 

As a last experiment, we tried Temperton's "best" scheme, in which heights 
were restored after each iteration during 10 iterations while velocities were 
allowed to adjust, followed by 10 iterati ‘s in which heights adjusted and winds 
were restored. We did obtain good results with this variation of Temperton's 
method, as shown by table 3. However this option still requires the equivalent 
of 120 0-R iterations, and fig. 7 shows clearly that in terras of balance, better 
results are obtained with just 12-15 0-R iterations. 

No experiments have been made with Meslnger's scheme, but the linear 
analysis in section 3 indicates that the results would be qualitatively similar 
to those of the NH2 method. 

b) Gradient wind approximation 

The main effect of the approximation of the wind by its geostrophic value 
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is to overestimate the velocity field In the regions with cyclonic curvature and 
uiiderestlmate It In the antlcyclonlc regions. Therefore we can make a simple 
analysis of the adjustment problem for a geostrophlc perturbation of the winds 
by assuming a perturbation of the form 

V ■ k X Vii I ill ■ asln k x sin ky, 4 “0. 

O /V o ’ ^o ’ ^o 

^'rom the linearized equation of conservation of potential vorticlty we can ob- 
tain the perturbation geopotential that remains in geostrophlc balance after 
inertia-gravity waves have been dissipated.: 

y2^ _ 1 ^ i - X (|. » 

from which 

K = ^ ^ ‘ “kRd)^^+ 1 

where Rj = — ^ is the Rossby radius of deformation. Computing the energy 
E * J/ (4> ^ + -I—) dx dy , we obtain 

E^/Eq - Ck-Rd)^/<^l^Rd)^ + 1} • 

This result indicates that in extratropical latitudes ^ where the scale 
of baroclinic waves is of the same order as the radius of deformation, most of 
the error introduced by the geostrophlc estimation of the wind will remain pre- 
sent after adjustment. This is also apparent in our numerical experiments: 
when tae heights were allowed to adjust, the reduction of rms wind error by the 
iterative schemes was small, while a sigi.if leant modification of the observed 
height field was introduced (see TABLE 3). The main discrepancy occured near 
the low, where the geostrophlc approximation strongly overestimated the winds 
and consequently the centrifugal force. During adjustment the height at the cen- 
ter of the low decreased from its original value of 2660 m to 2480 m, whereas 
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Che center of Che high Increased from 3150 m to 3160 m. 

We have seen Chat Che common practice of restoring Che mass fields after 
each Iteration Is not a solution to this problem, because It decreases enormously 
the efficiency of the dynamic Initialization schemes. A better Idea Is to decrease 
the Initial wind errors by the use of a more accurate approximation than the 
gcostrophlc for Che Initial winds. We found that a simple approximation of 
the gradient wind formula gave excellent results. 

The gradient wind equation can be written as 

V - V - ^ (5.2) 

g fr • 

where V Is the gradl^t wind, V Che geostrophlc wind and r the radius of 
curvature of the trajectory. The gradient wind approximation is strictly valid 
for circular steady flow but V approximates better than V the effects of flow 
curvature. However, the equation is quadratic in V and subject to a restriction, 
similar to the elllpticity condition, to ensure that the solution is real: 


+ > 0 . 

4 “ 


(5.3) 


If we assume that the right hand side of (5,2) Is small comp''red to V , 

S 

we can write V ■ V (1 + e) and neglect terms L_ e^. Then 

6 


e - - 


■A 


fr + 2V • 


(5.4) 


3/2 

The radius of curvature R = - (1 + y*) /y'* can be estimated from the 

geostrophlc streamlines; y’ is the slope of a curve Mx,y) = constant, and 

therefore y' ^ ^ , Combining these formulas we obtain 

X y 

/x2 . a2^3/2 

_ (tX Yy) ✓c c\ 

^ i *2 _ 2(j,% (j) + d, ‘ (5.5) 

^xx ^y ^x^y^xy ^xy^yy 


This formula allows the determination of r and V everywhere. In the 
regions where (5.3) is not satisfied, jel > 0.5, no longer a small number. 
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In these regions the correction of the geos trophic wind Is not Important because 

either V Is small > and so Is the error t or r Is large, and therefore curvature 
S 

effects are small. In our experiments we made no correction of the geostrophlc 
V i.:-d whenever |e| > ,5. 

This simple correction procedure was applied to the geostrophlcally per* 
turbed state. The resulting rms wind error was 3.8 m s ^ compared to 7.7 m s ' ^ 
In the geostrophlc error. When no further Initialization was made, gravity 
waves had an amplitude of 12 m compared to 125 m In the geostrophlc case. The 
0-R Iterative scheme was applied allowing heights to adjust, and again it con- 
verged In about 12 iterations , eliminating gravity waves . Table 3 indicated 

that the rms height error after Initialization wtas about 10 times smaller than 

5 

In the geostrophlcally perturbed case . However, these very encouraging results 
may be enhanced by the rather circular and steady state character of our basic 
flow (fig. 4). 

c) Initialization of the randomly perturbed state . 

In the previous experiments we Lave made no use of velocity measurements. 

In this section we describe a simulation of actual observations of both height 
and wind fields by perturbing the reference state fields with normally distributed 
random errors. The random errors are chosen to be compatible with real observa- 
tional errors. 

Three cases were run, with rms height errors of 0 m, 5 m and IG m respec- 
tively. In all cases an error of 3 m s ^ rms was introduced on each component 
of the velocity field. In the threa cases the gravity waves generated by the 
errors had an amplitude large enough, ri'f completely obliterate the reference fore- 
cast . 

Only the iterative schemes NH2 and 0-R were applied. The geopotential field 
was allowed to adjust and 150 complete iterations were performed. In both methods 
the rms error varied rapidly during the first 10 iterations and then attained 

5 

Lmilar results have been reported by Gauntlett and Seaman (1974). 
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a quasl-steavy value which continued to decrease very slowly as the Iterations 
proceeded. Again gravity wav«^s were completely eliminated by both methods. 

Table 4 presents the results and again we find that for the same n”’i.uer of 
complete iterations 0-R is somewhat faster than NH2, even though one NH2 itera- 
tion requires 3 times more computations than one 0-R Iteration. 

Another observation is that the magnitude of the height errors has no 
effect on the results. Again, this is explained by geostrophlc adjustment analy- 
sis. Impulsive perturbations on the height with a spaclal scale L will generate 
gravity waves which in an Infinite domain move away to infinity, and in our 
experiments are dissipated by the iterative schemes. The remaining geostrophlc 
mode contains most of the energy of the perturbation if L >> R^ » /$/f , and 
very little energy if L << R^ (see for example Chamey, 1973). In the case of 
observational errors their scale is the grid size, much smaller than the radius 
of deformation, and therefore most of the height perturbations goes into the 
dissipated gravity waves. The opposite situation occurs with velocity perturba- 
tions, and, as shown in table 4, a large part of the energy of wind errors 
remain in the fields after balance has been attained. In any case, the size of 
errors after initialization is reasonably small, 6 m for the heights and 2ms 
for the wind speed. 

With the balance equation approach only the heights were used. Case I, 
with no height errors, coincides with the case discussed in section 5a. In case 
II, we were able to "correct" the heights as described in section 5a, and satisfy 
the ellipticity condition everywhere, but the cycle-scan method of solution of 
the balance equation failed to converge. In case III the ellipticity condition 
was violated at many p< iu^s and the "correction" procedure failed to determine 
an elliptic geopotential field. Therefore we were unable to solve the balance 
equation in cases II and III. 

The ellipticity condition (5.1) can be written as 




rms ei 

height 

(m) 

rror 

wind^ 
(m sec 

1 

Method 

rms eri 
initii 
height 
(to) 

cor after 
allzation 
wind_ 

' (m sec”l) 

amplitude of 
gravity wave 
during 48 h 
forecast 

Case I 

0 

■ 

None 

— 

— 

40-350 

NH2 

6.4 

2.0 

0 

0-R 

6.2 

2.0 

0 

Balance 

0.1 

0.7 

3 

Case II 

5 

1 

None 

— 

— 

120-350 

NH2 

6.5 

1.9 

0 

0-R 

6.3 

1.8 

0 

Balance 

No convergenc 


Case III 

10 

1 

None 

— 

— 

70-350 

NH2 

6.6 

1.9 

0 

0-R 

6.5 

1.8 

0 

Balance i 

’’Correctiori* procedure failed 


Table 4. Random perturbations on the reference state 
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Ah « h - h < (5.6) 

In this £orm« It Is a restriction on the difference between the height at a 
point and the average height at the 4 adjacent points. It becomes increasingly 
stringent as the model resolution is Increased. For example, with As 250 km 
and f - 10~‘*sec Ah^^^ “ 8 ro, but if As = 125 km, the ellipticity condition 

is violated when Ah > 2 m! This is a. serious drawback of the balance equation 
approach, especially since the tendency of modem NWP models is towards smaller 
grid sizes. 
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6. Summary and conclusions 

We have presented a linear analysis and comparison of the damping properties 
of six Initialization schemes: the two proposed originally by Nltta and Hovermale 

(1969), the one by Meslnger (1972), a scheme attributed to Okamura (Nltta, 1969, 
Appendix), Temperton's (1973) averaging scheme and a modification of Okamura’ s 
scheme by Kalnay de Rivas. The linear analysis Indicates that the Okamura^Rlvas 
scheme has the most efficient damping properties over the whole frequency range, 
suggesting that It should be faster than the other methods and give more stable 
results . 

A nonlinear shallow water equations model on an f-plane has been used to 
test the initialization schemes and the results agree well with the linear 
analysis. When the heights are recovered after each Iteration, the iterative 
methods have a very slow convergence rate, because most of the Imbalarce is 
also recovered after each Iteration. When the heights are allowed to adjust 
freely, the iterative schemes converge much faster. In particular, the Okamura- 
Rivas scheme attains complete balance in only 12-15 iterations, equivalent to 
about 5 to 6 hours of regular forecasting using the leapfrog scheme. 

In the case of a geostrophic perturbation of the reference state, in 
which the observed winds are replaced by their geostrophic values, most of the 
error energy remains in the fields after free adjustment. This is in agreement 
with linear adjustment theory since the perturbation occurs in scales similar 
to the radius of deformation. 

Observational errors are also simulated by the introduction of random 
errors into the reference heights and velocity fields. In this case the dynamic 
initialization methods converge to a state much closer to the reference state 
than in the case of a geostrophic perturbation. In accordance with adjustment 
theory, small scale height errors are dissipated into gravity waves, while a 
significant portion of the small scale velocity error energy is retained in the 
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final fields. 

We conclude that whenever available, reliable wind observations should 
be Included In the Initial data for the dynamic initialization methods. The 
smoothing Incorporated In conventional analysis techniques may also help to 
reduce the errors In the observed fields. In data-s parse regions we suggest 
the use of a simple gradient wind approximation which can be directly evaluated 
from the geopotential field. In our experiments this approximation produced much 
better results than the geostrophlc wind approximation. Another method to 
Improve the Initial estimation of the wind field can be to solve the balance 
equation on a coarse grid In data-sparse regions and then Interpolate to finer 
resolution. In any case the best procedure is to first obtain a good estimate 
of the initial fields and then apply the iterative technique, allowing the free 
adjustment of the mass field. 

The balance equation approach provides good but not complete balance in 
the initial state for a primitive equation model. However this approach depends 
very critically on the ellipticlty condition, which in its simplest form is 
a restriction on the maximum amoxmt by which the height at a point can exceed 
the average height at the neighboring points. Uc have reported that this re- 
stiictLon is severely violated around strong anticyclones and, when the resolution is 
inereasedjby measurement errors typical of those occurring in atmospheric 
observations. The ellipticlty condition does not affect the Iterative techniques. 
When free adjustment is allowed the dynamic initialization methods become not 
only simpler but faster the balance equation approach.. In our numerical 

model, the Okamura-Rivas scheme required for convergence at least an order of 
magnitude less computations than the balance equation method. 

We have made our tests in a simple shallow water equations model, v;ith 
high frequency inertia-gravity waves. In baroclinic primitive equations models, 
inertia-gravity waves can have frequencies with values as low as the Coriolis 
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parameter. It Is In these models that the flexibility of the Okamura-Rlvas 
method may be put to maximum use. Even with the sequence n = 1, 1.6, 4, which 
is not optimum, and a time step of 15 minutes, waves with a frequency w = f = 

10 **8 ^ will be reduced by 20% in. 12 Okamura-Rivas iterations, and only by 

4% using the Nltta-Hovermale schemes and the same number of computations. 

The problem of four-dimensional data assimilation has not been considered 
here, but our results suggest that the use of the Okamura-Rivas scheme may be 
extremely useful in damping gravity waves generated by the introduction of new 
data in a model. A few Okamura-Rivas iterations after each set of data is 
Introduced followed, if necessary, by the use of a dissipative time scheme, will 
probably be sufficient to ensure a successful assimilat\on. 
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Figure Captions 

Figure 1 : Reduction of amplitude after one complete iteration for the NHl, 

NH2, and 0-R scheme. 

Figure 2 : Comparison of the relative efficiency of the NHl, NH2, 0 and 0-R itera- 

tive schemes. Reduction of amplitude after 12 F computations. 

Figure 3 ^ Same as Fig. 2 except for the Temperton, Mesinger and Okamura-Rlvas 
schemes. 

Figure 4 ; Reference height field. P indicates the point at which the height 
was monitored during the forecast. 

Figure 5 ; Reduction of error during initlaliKatlon. a) rms velocity error 

when the geopotential is restored after each iteration, b) rms velocity 
error during free adjustment, c) rms height error during free adjustment. 

Figure 6 ; Height forecast at the point P. Lower curves: heights restored 

after each iteration. Upper curves: heights adjusted freely. The 

reference forecast is also indicated. The initialization scheme used 
was a) Nitta-Hovermale 1, b) Nitta-Hovermale 2, c) Okamura, d) Okamura- 
Rivas . 

Figure 7 : Comparison of the height forecast at point P obtained after 12 and 

15 Okamura-Rivas iterations and after 20 "best" Temperton iterations. 
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